//////////////////////////////////////////////////////////////////////
//                                                                  //
//                                                                  //
//              Copyright (C), 2000  SEIKO EPSON Corp.              //
//              ALL RIGHTS RESERVED                                 //
//                                                                  //
//                                                                  //
//   file name : sinh.c                                             //
//                                                                  //
//   Revision history                                               //
//       00/02/10    first release                 M.Igarashi       //
//                                                                  //
//////////////////////////////////////////////////////////////////////

#include <math.h>
#include <smcvals.h>

#define	K0   0.336083852466226007246e0
#define	K1   0.3181332190110303793958e2
#define	K2   0.9178748106756360272984e3
#define	K3   0.6282813250833546179346e4
#define	S0  -0.1292607311299731203424e3
#define	S1   0.6282813250833546241924e4


// ALGORITHM
// 1.|x| > 0.5
//     sinh x = ( e^x - e^-x ) / 2
// 2.|x| <= 0.5 ( avoid the cancellation )
//     sinh x = x + x^3 / 3! + x^5 / 5!
//

//  the architecture of double floating point
//
//   0 1          11                  31 32                               63 bit
//   -----------------------------------------------------------------------
//  | |   exponent  |                    fraction                           |
//   -----------------------------------------------------------------------
//
//  |               |                   |                                   |
//  |    12 bits          20 bits       |             32 bits               |
//  |            lower word             |            higher word            |
//
//         bit    0         sign bit         (  1 bit  )
//              1 - 11      exponent part    ( 11 bits )
//             12 - 63      fraction part    ( 52 bits )
//
//

double sinh(double dfX){

	long lX;
	unsigned long ulLx,ulHx;
	double dfTemp,dfTemp2,dfX2,dfX3,dfRet;
	
	GETW_L(lX,dfX);		// get low
	GETW_H(ulHx,dfX);		// get high
	CLR_SIGN( dfX2, dfX );	// clear sign
	
	ulLx = lX&0x7fffffff;		// mask sign

	
	// |x| > 0.5
	
	if( ulLx > 0x40428000 || (ulLx == 0x40428000 && ulHx > 0x0 )){	// check qarea
		dfTemp = exp( dfX2 + LS );
		if ( (lX&0x80000000) == 0x0 ){			// positive
			dfRet = dfTemp;
			return dfRet;
		}else{									// negative
			dfRet = -dfTemp;
			return dfRet;
		}
	}
	
	if( ulLx > 0x3fe00000 || (ulLx == 0x3fe00000 && ulHx > 0x0 ) ){		//check area
		if ( (lX&0x80000000) == 0x0 ){			// positive
			dfRet=exp(dfX2+LS)-exp(-dfX2+LS);
			return dfRet;
		}else{									// negative
			dfRet=-(exp(dfX2+LS)-exp(-dfX2+LS) );
			return dfRet;
		}
	}
	
	// |dfX| <= 0.5
	
	dfX3=dfX2*dfX2;	
	dfTemp=K0*dfX3+K1;
	dfTemp=dfTemp*dfX3+K2;
	dfTemp=dfTemp*dfX3+K3;
	dfTemp*=dfX2;
	dfTemp2=dfX3+S0;
	dfTemp2=dfTemp2*dfX3+S1;
	dfTemp/=dfTemp2;

	if ((lX&0x80000000)==0x0){			// positive
		dfRet = dfTemp;
		return dfRet;
	}else{									// negative
		dfRet = -dfTemp;
		return dfRet;
	}

}
